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Octahedral Fe 2+ molecules are particularly interesting as they often exhibit a spin-crossover tran- 
sition. In spite of the many efforts aimed at assessing the performances of density functional theory 
for such systems, an exchange-correlation functional able to account accurately for the energetic 
of the various possible spin-states has not been identified yet. Here we critically discuss the issues 
related to the theoretical description of this class of molecules from first principles. In particular 
we present a comparison between different density functionals for four ions, namely [Fe(H20)e] 2+ , 
[Fe(NH 3 ) 6 ] 2+ , [Fe(NCH) 6 ] 2+ and [Fe(CO) 6 ] 2+ . These are characterized by different ligand-field 
splittings and ground state spin multiplicities. Since no experimental data are available for the gas 
phase, the density functional theory results are benchmarked against those obtained with diffusion 
Monte Carlo, one of the most accurate methods available to compute ground state total energies 
of quantum systems. On the one hand, we show that most of the functionals considered provide 
a good description of the geometry and of the shape of the potential energy surfaces. On the 
other hand, the same functionals fail badly in predicting the energy differences between the various 
spin states. In the case of [Fe(H 2 0) 6 ] 2+ , [Fe(NH 3 ) 6 ] 2+ , [Fe(NCH) 6 ] 2+ , this failure is related to the 
drastic underestimation of the exchange energy. Therefore quite accurate results can be achieved 
with hybrid functionals including about 50% of Hartree-Fock exchange. In contrast, in the case of 
[Fe(CO)e] 2+ , the failure is likely to be caused by the multiconfigurational character of the ground 
state wave-function and no suitable exchange and correlation functional has been identified. 



INTRODUCTION 

Among the transition metal complexes, octahedral 3d 6 
Fe 2+ molecules are systems of particular interest. In fact 
they often undergo the spin-crossover (SC) transition [1] . 
In their most common form the molecules ground state is 
a spin singlet (S — 0), with the Fe six 3d electrons paired 
up in the t2 g ir* antibonding orbitals. Upon increasing 
the temperature, the high spin quintet state (5 = 2), in 
which two electrons are promoted from the tig^* to the 
e s cr* orbitals, becomes thermodynamically more stable 
(see Fig. [T] for the molecular orbital diagram) . The SC 
transition is entropy driven and it is regulated by the 
difference in the Gibbs free energy of the HS and LS 
states, 

AG = G HS - G LS = AH - TAS . (1) 

Here AH = i?ns — #ls an d AS — Shs — Sls indicate 
respectively the enthalpy and the entropy variation (note 
that for AG > the LS configuration is more thermo- 
dynamically stable than the HS one). For SC molecules 
AH > 0, but also S H s > S L s, i.e. AS > 0. Hence, 
for large enough temperatures (T > T c = AH/ AS), the 
cntropic term dominates over the enthalpic one and the 
molecules transit from the LS to the HS configuration. 
There are two contributions to the entropy: the first is 
provided by the spin and the second by the molecule vi- 
brations. In fact, the transfer of two electrons to the e g a* 
orbitals, which are more-antibonding than the tigit* , 
weakens the chemical bond and produces a breathing of 



the metal ion coordination sphere. This results in a soft- 
ening of the phonon modes and then an increase of the 
vibronic entropy. 

The SC transition is usually reported either for 
molecules in solution or in single crystals and, depending 
on the strength and on the origin of the inter-molecular 
interactions, it can have cooperative nature and present a 
thermal hysteresis loop. Interestingly, the transition can 
be also induced by illumination. This phenomenon is 
called LIESST effect (Light-Induced-Excited-Spin-State- 
Trapping) and it is explained through a mechanism in- 
volving several excited states [2]. The lifetime of the 
metastable HS state is usually very long at low tem- 
perature as the relaxation to the ground state is due 
to the weak electronic coupling between these states [3]. 
The opposite process, in which molecules populating the 
metastable HS state are brought back to their ground 
state, is also possible and it is called reverse LIESST ef- 
fect. 

Although, SC molecules have been traditionally stud- 
ied for possible applications in storage and sensor devices 
[H El H] i they have recently emerged as promising materi- 
als for molecular spintronics [6-8]. In fact, the electronic 
transport through these molecules has been predicted 
[HI HO] and experimentally reported [TTrTB] to depend 
strongly on their spin state. 

Given such renewed interest in spin crossover com- 
pounds there is also a growing fundamental effort in mod- 
eling their properties. In this respect one aims at using 
an electronic structure theory, which is at the same time 
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FIG. 1: Energy level diagram for an octahedrally coordinated 
transition metal (TM) ion. In crystal field theory, the 3d or- 
bitals of the TM ion have an energy split Acf due to the 
electrostatic interaction with the ligands. In contrast, in lig- 
and field theory the 3d orbitals of the TM ion form covalent 
bonds with the ligands. In this diagram we assume that each 
ligand contributes three p-orbitals, one with the positive lobe 
pointing toward the TM ion, cr-type, and two with the lobes 
perpendicular to it, 7r-type, (note that the a- and 7r-type p 
orbitals are degenerate but here they are plotted slightly sep- 
arated for better display). The 7r-type p orbitals couple with 
the TM t2 g states, while the cr-type p orbitals couple with 
the e g . Since the ir interaction is weaker than the a one, the 
antibonding t2 g ir* orbitals lie lower in energy than the e 3 cr* 
ones. The energy splitting between the £2 9 7r* and the e g a* 
orbitals is indicated by Alf + Acf. 



accurate and scalable. Accuracy is needed for reliable 
predictions of the spin crossover temperature, while scal- 
ability is required by the size of the typical molecules. 
This becomes particularly crucial for molecules in in crys- 
tals and when deposited on metallic surfaces, since the 
typical simulation cells are large. Density functional the- 
ory (DFT) is in principle both scalable and accurate, but 
to date it is completely unclear how it does perform rel- 
atively to this problem. 

In this paper we investigate the performance of DFT 
against four model Fe 2+ ions, namely [Fe(H20)e] 2+ , 
[Fe(NH 3 ) 6 ] 2 +, [Fe(NCH) 6 ] 2 + and [Fe(CO) 6 ] 2+ . In par- 
ticular we answer to the question on which exchange- 
correlation functional performs best. As no experimental 
data are available for these ions we benchmark against 
diffusion Quantum Monte Carlo (QMC) results. Note 
that often experimental data are difficult to compare with 
microscopic theory, so that our work provides a rigorous 
benchmark for the theory itself, which is probably more 
informative than a simple comparison to experiments. 
The paper is organized as follows. In the next section 
we present an overview of the problem and of the var- 
ious known shortcomings of DFT, and we discuss criti- 
cally which elements one has to consider when comparing 
electronic structure data to experiments. Then we will 
provide some computational details and move to the re- 



sults. First we will discuss our DFT calculations for the 
four different ions and then we will compare them with 
the QMC ones. Finally we will conclude. 




FIG. 2: Potential energy surface of the HS and LS state of 
a SC molecule. The collective coordinate r represents all 
the 3N nuclear coordinates of the molecule. The zero point 
phonon energies for the HS and LS state, -Ehs E and £ls E , 



the adiabatic energy gap 
gaps 

indicated 



AEIT = A£; vcrt (r L s) 



A_B a , and the vertical energy 
and A-Egjf = AE vert (r HS ) are 



DENSITY FUNCTIONAL THEORY FOR 
SPIN-CROSSOVER MOLECULES: STATE OF 
ART 

When we consider a single molecule in vacuum at zero 
temperature AG coincides with the internal energy dif- 
ference, which reads in the adiabatic approximation as 



AE = AE 



adi; 



AE 



ZPE 



(2) 



Here, AE ZPE = E^ E - E^ E and E^ E (£ff E ) is the 
zero-point phonon energy of the HS (LS) state, while 



A £ad Ia = £ Hs(rHg) _£ Ls(rLs) 



(3) 



is the adiabatic energy difference (often called "adiabatic 
energy gap"). The collective coordinate r represents the 
3N nuclear coordinates of the molecule and the energies 
^Hs( r ) an d ^Ls( r ) define the potential energy surfaces 
(PESs) respectively of the HS and LS state (see Fig. |2]). 
In addition to the adiabatic energy gap we can also define 
the vertical energy difference ("vertical energy gap" ) [T5] 



AE vert (r)=E m (r)-E hS (r) 



(4) 



and the difference of vertical energy gaps (DOG) 16J: 

DOG(r 2 ,n) = A£ vcrt (r 2 ) - AE VCIt ( n ) . (5) 

All of these quantities can be computed by using ab-initio 
electronic structure methods. As we mentioned in the 
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introduction, DFT is always the preferred one. In fact 
SC molecules are composed of, at least, fifty atoms and 
a good balance between expected accuracy and compu- 
tational cost is required. However, there are many is- 
sues connected to the DFT description of SC molecules, 
which either have not been properly addressed or have 
not found any satisfactory solution yet. Here we list some 
of them. 

- The Functional dilemma. For each Fe 2+ 
molecules, in general, every exchange-correlation 
functional returns a very different adiabatic energy 
gap (see for example references (TBI HI])- These 
differences can be as large as few eV and different 
functionals do not sometimes even predict the same 
A£' adla sign. In a nutshell, no agreement around 
which functional performs best has been reached 
so far (the discussion below will explain how prob- 
lematic is a direct assessment of the DFT results 
through a comparison with the experimental data). 
However, some general trends, which relate func- 
tionals belonging to the same "class" (or the same 
"rung" if we refer to the "Jacob's ladder" [TB] clas- 
sification scheme), have been identified. 

1) the local density approximation (LDA) (first 
rung) tends to stabilize the LS state. This is due to 
the underestimation of the exchange energy |19j . In 
particular, the exchange keeps electrons of like-spin 
apart so that their Coulomb repulsion is reduced. 
Therefore, the exchange underestimation is accom- 
panied by the overestimation of the Coulomb en- 
ergy for two electrons of equal spin. This, in turns, 
leads to the stabilization of the LS state at the ex- 
pense of the HS state. 

2) the generalized gradient approximation ( GGA ) 
and the meta-GGAs (second and third rungs) \2U§ 
give results that depend on the specific compound 
and on the exact DFT conditions that each func- 
tional satisfies. "Traditional" GGAs, such as 
the Perdew-Burke-Ernzerhof (PBE) [5T] and the 
BLYP, which combines the Becke exchange [22] 
with the Lee- Yang-Parr (LYP) correlation [23], re- 
duce only slightly the LDA over-stabilization of 
the LS state. Therefore they do not represent a 
drastic improvement. In contrast, some more re- 
cent GGA functionals, such as the OLYP (com- 
bining the Handy and Cohen's OPTX exchange 
[2"1] with the LYP correlation), have been claimed 
to perform rather well [55] [5^1 [25] ■ Among the 
meta-GGAs, the Van Voorhis-Scuseria exchange- 
correlation (VSXC) functional [57] was tested by 
Ganzenmuller et al. 28], who concluded that 
it provided quite accurate results for single-point 
calculations, while it predicted artificially twisted 
structures. 

3) hybrid functionals (forth rung) tend to favor the 



HS with respect to LS one. Reiher and co-workers 
[29l [30] argued that the amount of Hartree-Fock 
(HF) exchange in many hybrid functionals is too 
large to predict correct energy gaps. They then 
proposed a re-parametrization of the B3LYP func- 
tional 3T], called B3LYP*, which includes only 15% 
of HF exchange (in contrast to the standard B3LYP 
20%). Although B3LYP* is believed to give satis- 
factory results for several SC complexes, some stud- 
ies suggested that a further reduction of the HF ex- 
change could be needed in order to describe others 
Fe 2+ compounds 28, Mj- hi contrast, the amount 
of HF exchange in B3LYP was judged insufficient 
for the small ions [Fe(H 2 0) 6 ] 2+ and [Fe(NH 3 ) 6 ] 2 +. 
For these systems it has been claimed that PBEO 
[3HE3]j which includes up to 25% of HF exchange, 
gives more satisfactory results 25, 26J. In practice, 
for each compound, the amount of HF exchange 
can be varied to fit the desired values for the gaps 
but no "universally good choice" has been identi- 
fied so far. The dependence of the adiabatic en- 
ergy gap on the amount of HF and either of LDA 
(Slater) or GGA (B88 and OPTX) exchange is ex- 
plained very clearly in the work of Ganzenmuller 
et al. 28 . Finally it is important to remark that, 
even when a certain hybrid functional is found to 
return satisfactory energy gaps, it might not be the 
optimal functional to describe other properties of 
the molecule. 

- Basis set. DFT calculations for SC molecules 
are usually performed by using Quantum Chem- 
istry codes, where the wave-function is expanded 
over either Gaussian (GTOs) or Slater-type orbitals 
(STOs) p2]. In many cases, the values of the en- 
ergy gaps depend substantially on the choice of the 
basis sets and on the spatial extension of the local 
orbitals [35]. Although this is a less severe prob- 
lem compared to that of identifying the exchange- 
correlation functional, it must be kept into consid- 
eration. In principle, the use of plane-waves basis 
sets [36], instead of GTOs and STOs, could be a so- 
lution, but in practice plane-wave calculations are 
are computationally expensive because of the need 
of satisfying periodic boundary conditions. 
A large number of plane-waves is usually needed 
as the electronic density is concentrated in a small 
fraction of the supercell volume. Furthermore very 
large supercells are typically required. This is due 
to the fact that SC complexes are often 2+ ions. 
Therefore a negatively charged background must be 
introduced in order to preserve the overall charge 
neutrality so that the total energy remains bound. 
The energy calculated in this way approaches then 
the one for an isolated system only in the limit 
of large supercell and, unfortunately, such conver- 
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gence is slow. Although corrections to the expres- 
sion of the computed energy have been proposed 
[37] [38] , this effect can be properly accounted for 
only by considering large supercells and by per- 
forming multiple calculations for supercells of dif- 
ferent sizes. 

- Geometry optimization. Each class (rung) of 
exchange-correlation functionals returns different 
metal-ligands bond-lengths. Usually, LDA gives 
shorter bonds than hybrids functionals, while stan- 
dard GGA results are in between these two ex- 
tremes E0]- Although these differences are 
usually quite small, less then 0.1 A against an 
average bond-length of about 2 A, they might 
strongly affect the electronic properties. Indeed a 
very delicate balance between ligand-field splitting 
and Hund's coupling establishes the spin state (see 
Fig. [I). Then, small errors in the predicted geome- 
tries can drastically alter this balance. 
Unfortunately, the quality of the DFT relaxation 
can not be easily assessed through a comparison 
with available experiments. In fact, while DFT 
calculations are usually carried out for molecules 
in the gas phase, the experimental geometries are 
obtained through X-ray measurements for crystals 
|41) . As far as we know, no detailed DFT studies 
about SC molecules in the condensed phase have 
been published so far. Furthermore, such a study 
must face the additional difficulty of the need of 
accounting for inter-molecular interactions. These 
have usually dispersive nature and, therefore, they 
are either not described or badly-described by most 
of the popular functionals [42]. 

- Electrostatic contributions to the total en- 
ergy. SC molecules, which are usually in the 2+ 
charging state, are surrounded by counter-ions (for 
example [PFg] 2 ~ or [BF<4] 2 ~). Because of the elec- 
trostatic potential generated by such counter-ions 
and by the other molecules in the crystal, the to- 
tal energy of a SC complex in the condensed phase 
differs from that in the gas phase. Furthermore, 
since at the SC phase transition there is a charge 
redistribution over the individual molecules and a 
lattice expansion, such electrostatic potential does 
not induce a simple rigid shift of the minima of the 
HS and LS PESs (compared to those calculated for 
the gas phase). The energy gap then turns out 
to be different for the same molecule in different 
phases [HI [H] . Unfortunately this effect is always 
neglected by DFT calculations, which aim at as- 
sessing the performances of DFT by using experi- 
mental data. 

- Finite-temperature effects. Special care should 
be taken in order to include finite-temperature ef- 



fects when comparing DFT to experiments. In fact, 
at finite temperature instead of the adiabatic en- 
ergy gap, the Gibbs free energy, Eq. |l]), must be 
computed. So far, this has been attempted only by 
Ganzenmuller et al. |28j . However, unfortunately, 
their calculations did not fully account for either 
the electrostatic contribution to the total energy or 
the effect of the periodic lattice on the molecular 
structure. 

This list clearly emphasizes how the main handicap 
in the theoretical description of SC complexes is related 
to the difficulties in assessing the performances of any 
given exchange-correlation functional. Any benchmark 
involving a comparison against experimental results is 
fated to fail, unless vibrational, environmental, crystal- 
lographic and finite-temperature effects are properly ac- 
counted for. However, this task is at present too demand- 
ing to be practically achievable. In contrast, as pointed 
out by Fouqueau et al. [25] [26] , the current best strategy 
consists in providing benchmark values for various inter- 
esting quantities through highly accurate ab-initio meth- 
ods. These can then be compared with the DFT results 
in order to identify which functional performs better. 

In some interesting works [21 [25] [23 HSJ US], wave- 
function based methods were considered (see discussion 
below). However, unfortunately, the authors themselves 
admitted that their results were plagued by a number of 
systematic errors. These were ascribed to the too small 
basis set used for the Fe 2+ ion and to the fact that such 
computational methods describe exactly only for static 
electronic correlation, but do not include the dynamic 
one, which, however, can be partially included by per- 
turbation theory. 

Here we have chosen to employ as benchmark electronic 
structure theory diffusion Monte Carlo (DMC) [47H49] . 
This represents one of the most accurate computational 
techniques to calculate the energy of a quantum system 
and it is able to return up to 98% of the correlation en- 
ergy (including both static and dynamic contributions). 

In this work, we compare systematically several DMC 
results to those obtained with DFT for a few selected 
Fe 2+ complexes. Unfortunately, such a systematic in- 
vestigation requires a large use of computational re- 
sources and, therefore, it can not be carried out for 
molecules composed of tens of atoms (such as the most 
typical SC complexes). We have then focused our atten- 
tion on the following ions: [Fe(H 2 0) 6 ] 2+ , [Fe(NH 3 ) 6 ] 2+ , 
[Fe(NCH) 6 ] 2 + and [Fe(CO) 6 ] 2+ , which are shown in Fig. 
[3] These are small enough to allow several DMC calcu- 
lations to be performed at a reasonable computational 
cost. Furthermore, and more importantly, the study of 
their electronic structure presents all the problems men- 
tioned above so that our main conclusions can be ex- 
tended to large SC molecules as well. Finally, according 
to the spectrochemical series [SU] , these ions have a differ- 
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ent ligand field splitting. The one of [Fe(H 2 0)e] 2+ is the 
smallest and that of [Fe(CO)g] 2+ the largest. We then ex- 
pect that Fe(H 2 0) 6 ] 2 +, [Fe(NH 3 ) 6 ] 2+ and [Fe(NCH) 6 ] 2+ 
have a HS ground state, while [Fe(CO) 6 ] 2+ a LS one. 
Therefore, our study scans through systems of different 
ground state spin multiplicity, it reveals general trends 
and it points to the systematic errors of each DFT func- 
tional. 



COMPUTATIONAL DETAILS 

DFT calculations are performed with the NWCHEM 
package [5T]- We use several functionals belonging to 
different "classes": 1) the default LDA of Vosko-Wilk- 
Nussair [5J], 2) the GGA BP86 functional, which com- 
bines the Becke88 exchange functional [22] with the Pre- 
dew86 correlation one [53] and 3) the hybrid functionals 
B3LYP [31], PBEO [31133] and the Becke "half and half" 
(Becke-HH) [53]. These include respectively 20%, 25% 
and 50% of HF exchange (note that, in the NWCHEM 
package, only an approximate version of the HH-Becke 
functional is currently implemented [55]). 

We have chosen three different basis sets: 1) 6-31G* 
(basis set called A), 2) the LANL2DZ basis set and pseu- 
dopotential [55] for Fe combined with the basis set 6- 
31++G** for all other atoms (basis set B) and 3) the 
triple-zeta polarized basis set of Ahlrichs [57] (basis set 
C). Geometry optimizations are performed both with and 
without specifying the molecule point group. The geome- 
tries optimized with these two strategies usually return 
consistent results with bond-lengths differences, which 
are within ±0.005A. We always check that the phonon 
frequencies are all real so that the final geometries corre- 
spond to stable energy minima. 

DMC calculations are performed by using the CASINO 
code [48]. The imaginary time evolution of the 
Schrodinger equation has been performed with the usual 
short time approximation and time-steps of various sizes 
are considered (typically At = 0.0125,0.005,0.001 a.u.). 
Expect for few cases, the energy differences are usually 
found to be converged with respect to the time-step errors 
already for At — 0.0125 a.u. Calculations are performed 
by using Dirac-Fock pseudopotentials [55] [59] with the 
"potential localization approximation" (PLA) [BB]. For 
[Fe(NCH)g] 2+ DMC simulations with this approximation 
are found to be unstable as the number of walkers "ex- 
plodes". Therefore we have used instead the "T-move" 
scheme [BT] [55], which eliminates the need of the PLA 
and treats the non-local part of the pseudopotential in a 
way consistent with the variational principle. The sim- 
ulations then become more stable. The single-particle 
orbitals of the trial wave function are obtained through 
(LDA) DFT calculations performed with the plane-wave 
(PW) code QUANTUM ESPRESSO [BJ. The same 
pseudopotentials used for the DMC calculations are em- 



ployed. The PW cutoff is fixed at 300 Ry and the PW 
are re-expanded in terms of B-splines [64] . The B-spline 
grid spacing is a = 7r/G max , where G max is the length 
of the largest vector employed in the PW calculations. 
Periodic boundary conditions are employed for the PW- 
DFT calculations and supercells as large as 40 A are 
considered. In contrast, no periodic boundary conditions 
are imposed for the DMC simulations. 

DMC calculations are performed for the molecular ge- 
ometries previously optimized by DFT. Therefore we can 
compare the DMC energies of molecular structures ob- 
tained by employing different functionals and basis sets. 
However, for each system, these energy differences are 
often smaller than the computed error bars. Only LDA 
systematically returns molecular geometries with much 
higher DMC energies than those obtained by using ei- 
ther GGA or hybrid functionals. 

DFT RESULTS 

[Fe(H 2 0) 6 ] 2+ 

The lower energy geometry of [Fe(H 2 0)g] 2+ is found 
to have C, symmetry for both BP86 and hybrid func- 
tionals. This is consistent with the results of Pierloot et 
al. [B5J. In contrast, with LDA, we were able to obtain 
relaxed atomic positions for both the HS and LS states 
only by using the the basis set A and without specifying 
the molecule point group. 

As expected from our introductory discussion, the 
molecule in the LS state has metal-ligand bond-lengths 
shorter than those of the molecule in HS state (by about 
7%). However, the details of the geometry depend on 
both the functional and the basis set. This can be clearly 
seen by inspecting Tab. [T] which reports a full list of the 
calculated Fe-0 bond-lengths for both the HS and LS 
states. On the one hand, LDA overbinds the molecule as 
compared to GGA and hybrids. On the other hand, the 
basis set A tends to shrink the Fe-0 bond-lengths when 
compared to the basis set B and C. Although the choice of 
basis set does not affect the bond-lengths as drastically 
as the functional does, it still influences the geometry 
greatly. The calculations performed with the basis set A 
return a quite large inclination (about 5 degrees) of the 
O-Fe-0 axis with respect to the 90 degrees angle it forms 
with the equatorial plane of the molecule. Furthermore, 
for basis set A and B, either the axial waters, which form 
the ligands, move "in" and the equatorial waters move 
"out" of their plane or viceversa. These results do not 
depend qualitatively on the functional. In contrast, the 
inclination of the O-Fe-0 axis and the distortion of axial 
waters disappears when the calculations are carried out 
by using the basis set C. 

Tab. [IT] shows our calculated values for the adiabatic 
energy gaps, where a positive (negative) energy means 
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FIG. 3: (Color on line) The cations investigated in this work [Fe(H 2 0) 6 ] 2+ (a), [Fe(NH 3 ) 6 ] 2+ (b), [Fe(NCH) 6 ] 2+ (c) and 
[Fe(CO)e] 2+ (d). Color code: C=yellow, 0=red (small sphere), Feared (large sphere), N=grey, H=blue. 



that the LS state has lower (higher) energy than that of 
the HS one. For each functional, our results are always 
in very close agreement with those obtained by other au- 
thors [SS] (the results are presented in cm -1 as well 
as in eV in order to allow for a better comparison with 
the various values found in literature). Here we can dis- 
tinguish a clear trend, summarized by the series: 

-A£ adia (LDA) < -A£ adia (GGA) < -A£ adia (B3LYP) 
< -A.E adia (PBE0) < -A£ adia (HH) . (6) 

This suggests that the calculated Ai? adla is strictly re- 
lated to the amount of HF exchange incorporated into 
the given functional. By increasing such contribution, 
we systematically stabilize the HS configuration with re- 
spect to the LS one. 



Functional Basis set 


A£ adia (cm" 1 ) 


A£ ,adia ( e y) 


LDA 


A 


-3986 


-0.4942 


BP86 


A 


-8989 


-1.1145 


BP86 


B 


-8381 


-1.0391 


BP86 


C 


-8400 


-1.0415 


B3LYP 


A 


-11589 


-1.4369 


B3LYP 


B 


-11027 


-1.3672 


B3LYP 


C 


-11045 


-1.3694 


PBEO 


A 


-14670 


-1.8189 


PBEO 


B 


-15512 


-1.9233 


PBEO 


C 


-14045 


-1.7414 


HH 


C 


-19620 


-2.4326 


HH 


C 


-18223 


-2.2594 



TABLE II: Adiabatic energy gap, A£? adia , for the cation 
[Fe(H20)e] 2+ - The functional and the basis set used for the 
each calculation are indicated. 



Functional Basis set 


rf LS (A) 


dm (A) 


LDA 


A 


1.917 


2.052,2.083,2.057 


BP86 


A 


1.985 


2.152,2.151,2.111 


BP86 


B 


2.02 


2.174,2.164,2.132 


BP86 


C 


2.01 


2.161,2.155,2.125 


B3LYP 


A 


2.005 


2.152,2.152,2.111 


B3LYP 


B 


2.003 


2.146,2.157,2.112 


B3LYP 


C 


2.029 


2.172,2.16,2.137 


PBEO 


A 


1.99 


2.1,2.145,2.134 


PBEO 


B 


2.013 


2.168,2.156,2.129 


PBEO 


C 


2.008 


2.152,2.147,2.124 


HH 


B 


2.010 


2.168,2.133,2.132 


HH 


C 


2.008 


2.149,2.131,2.128 



TABLE I: Bond-lengths of [Fe(H 2 0) 6 ] 2+ in the HS and LS 
state, as calculated with various functionals and basis sets. 
Note that the LDA calculations for the HS state did not cov- 
erge in the case of the basis set B and C. 



[Fe(NH 3 ) ( 



12+ 



The optimized structure of the [Fe(NH3)g] 2+ ion, cal- 
culated either with BP86 or with hybrid functionals, has 
a Z?3 symmetry for the LS state. This is further low- 
ered to C2 for the HS one. Our results are again consis- 
tent with those of Pierloot et al. [55] ■ Like in the case 
of [Fe(H 2 0)6] 2+ , we were not able to find the relaxed 
atomic geometry with LDA. Even when the geometry 
optimization procedure converges, like in the case of the 
basis set A and C, the minimum is found to be unsta- 
ble. This is indicated by the negative eigenvalues of some 
phonon modes. Nevertheless we report these results for 
completeness. 

The molecule in LS state has shorter average Fe- 
ligand bon d-lengths than the molecules in HS state (see 
III). In contrast to the case of [Fe(H 2 0) 6 ] 2+ , 
l3)e] 2+ does not show any strong deviation of the 



Tab 
[Fe(NH 

N-Fe-N axis with respect to the axis normal to the equa- 
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Functional Basis set 


d LS (A) 


dns (A) 


LDA 


A 


1.942 


2.188,2.162,2.160 


LDA 


C 


1.995 


2.204,2.201,2.214 


BP86 


A 


2.026 


2.267,2.254,2.253 


BP86 


B 


2.078 


2.30,2.295,2.274 


BP86 


C 


2.085 


2.279,2.302,2.289 


B3LYP 


A 


2.076 


2.281,2.281,2.275 


B3LYP 


B 


2.114 


2.315,2.296,2.294 


B3LYP 


C 


2.122 


2.32,2.308,2.283 


PBEO 


A 


2.05 


2.254,2.256,2.256 


PBEO 


B 


2.082 


2.292,2.277,2.272 


PBEO 


C 


2.093 


2.284,2.294,2.263 


HH 


c 


2.11 


2.296,2.286,2.266 



TABLE III: Bond-lengths of [Fe(NH 3 ) 6 ] 2+ in the HS and LS 
state, as calculated with various functionals and basis sets. 
Note that the LDA calculations for the HS state did not cov- 
erge in the case of basis set B. 



Functional Basis set c^ls (A) dns (A) 



LDA C 1.854 2.066,2.067,2.11 

BP86 C 1.917 2.171,2.171,2.155 

B3LYP C 1.974 2.206,2.201,2.201 

PBEO C 1.950 2.194,2.181,2.181 

HH C 1.990 2.20,2.196,2.196 



TABLE V: Bond-lengths of [Fe(NCH) 6 ] 2+ in the HS and LS 
state as calculated with various functionals and for the basis 
sets C. 



Functional Basis set A£ adia (cm" 1 ) A£ adia (eV) 



LDA 


C 


19126.3 


2.37135 


BP86 


C 


8410.89 


1.04282 


B3LYP 


C 


-1667.48 


-0.20674 


PBEO 


C 


-3544.58 


-0.43947 


HH 


C 


-12029.62 


-1.49148 



Functional Basis set 


A£ adia (cm" 1 ) 




LDA 


A 


8937 


1.1081 


LDA 


C 


7746 


0.9605 


BP86 


A 


195 


0.0242 


BP86 


B 


708 


0.0878 


BP86 


C 


672 


0.0834 


B3LYP 


A 


-5312 


-0.6586 


B3LYP 


B 


-4007 


-0.4969 


B3LYP 


C 


-4738 


-0.5874 


PBEO 


A 


-7695 


-0.9541 


PBEO 


B 


-7665 


-0.9504 


PBEO 


C 


-7117 


-0.8825 


HH 


c 


-13556 


-1.68077 



TABLE IV: Adiabatic energy gap, A£ adia , for the cation 
[Fe(NH3)g] 2+ . The functional and the basis set used for each 
calculation are indicated. 



torial plane for any combination of functionals and basis 
sets. 

Tab. |IV| shows several values for the adiabatic energy 
gaps. Once again these can be ordered according to the 
series Here, the LDA adiabatic energy gap indicates 
that [Fe(NH3)e] 2+ is LS. This result is even qualitatively 
incorrect as this cation is known to be stable in the HS 
state. BP86 also predicts the LS state to be the lowest 
in energy, although the value of AE adm is very small 
and probably very sensible to the exact details of the 
calculation. In fact, in contrast to our results, which 
agree with those in Ref. [26| > Pierloot et al. [65] obtained 
a negative value equal to about —0.2 eV. Finally hybrid 
functionals predict the ground state to be HS with the 
value of the gap being proportional to the amount of HF 
exchange included in the functional. 



TABLE VI: Adiabatic energy gap, A£ adia , for the cation 
[Fe(NCH) 6 ] 2+ calculated with various functionals and the ba- 
sis sets C. 



[Fe(NCH)e 



The results for the [Fe(H 2 0) 6 ] 2+ and [Fe(NH 3 ) 



12 + 



ions demonstrate that, although a good choice of basis set 
might be important to predict accurate molecular struc- 
tures, the estimated values of the adiabatic energy gaps 
depend mainly on the functional used. Indeed, for a given 
functional, two adiabatic gaps obtained with two differ- 
ent basis sets differ at most by a few tens of meV. This 
has to be compared with the differences in the values 
predicted by different functionals, which can be of sev- 
eral hundreds meV. For the ion [Fe(NCH)g] 2+ we have 
then decided to compare only calculations performed us- 
ing the basis set C, which typically gives us the lowest 
energy. 

[Fe(NCH) 6 ] 2+ has perfect octahedral symmetry in the 
LS state. In contrast, the structure of the HS state is pre- 
dicted to have D±h symmetry by B3LYP and PBEO and 



Ci symmetry by BP86 and Becke-HH. Tab. VI displays 
the values of adiabatic energy gap calculated with each 
functional. Once again these can be ordered according 
to the series (|6| . We find that the total energy of the LS 
state is at least 1 eV lower than that of the HS state for 
both the LDA and BP86. In contrast PBEO and B3LYP 
return the HS state as the most stable, but the abso- 
lute value of Ai? adla is only a few hundreds meV (note 
that our B3LYP adiabatic energy gap is consistent with 
that calculated by Bolvin Finally, the Becke-HH 

predicts A£ adia « -1.5 eV. 
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Functional Basis set 


c(ls (A) d H s (A) 


System 


Functional AE ZPE (cV) AE 


ZPE (cm" 1 ) 


LDA 


c 


1.848 2.199,2.172,2.123 


[Fe(H 2 0) 6 ] 2+ 


BP86 


—0 08195 


-661 


BP86 


C 


1.900 2.226,2.331 


[Fe(H 2 0) 6 ] 2+ 


B3LYP 


-0.09079 


-732 


B3LYP 


c 


1.948 2.307,2.367 


[Fe(H 2 0) 6 ] 2+ 


PBEO 


-0.09308 


-750 


PBEO 


c 


1.915 2.276,2.345 


[Fe(H 2 0) 6 ] 2+ 


HH 


-0.10272 




HH 


c 


1.915 2.322,2.329,2.366 


[Fe(NH 3 ) 6 ] 2+ 


BP86 


-0.17413 


-1404 








[Fe(NH 3 ) 6 ] 2+ 


B3LYP 


-0.16099 


-1298 


TABLE VII: Bond-lengths of [Fe(CO) 6 ] 2+ in the HS and LS 
state, as calculated with various functionals and for the basis 
set C. 


[Fe(NH 3 ) 6 ] 2+ 
[Fe(NH 3 ) 6 ] 2 + 
[Fe(NCH) 6 ] 2 + 


PBEO 
HH 

BP86 


—0 1 7636 
—0.16071 
—0 165Q3 


-1422 
-0.16071 
-1338 








[Fe(NCH) 6 ] 2+ 
[Fe(NCH) 6 ] 2+ 
[Fe(NCH) 6 ] 2 + 


B3LYP 
PBEO 


-0.13846 


-1239 
-1116 


Functional Basis set 


A£; adia (cm" 1 ) A£ adia (eV) 


LDA 


c 


41148 5.1017 


HH 


-0.14487 


BP86 


c 


27575 3.4189 


[Fe(CO) 6 ] 2 + 


BP86 


-0.20800 


-1677 


B3LYP 


c 


10656 1.32126 


[Fe(CO) 6 ] 2+ 


B3LYP 


-0.19894 


-1604 


PBEO 


c 


10888 1.3501 


[Fe(CO) 6 ] 2+ 


PBEO 


-0.21157 


-1706 


HH 


c 


5232 -0.6488 


[Fe(CO) 6 ] 2 + 


HH 


0.21799 





TABLE VIII: Adiabatic energy gap, A£ a , for the cation 
[Fe(CO)e] 2+ calculated with various functionals and the basis 
set C. 



[Fe(CO) e 



The [Fc(CO)g] 2+ ion has perfect octahedral symmetry 
in the LS state. This is then reduced to D±h in the HS 



(the metal-ligand bond-lengths are listed in Tab. VII I. 
The calculated adiabatic energy gaps are displayed in 



Tab. |VIII[ ). Again LDA and BP86 are found to (mas- 
sively) over-stabilize the LS state and the adiabatic en- 
ergy gap turns out unrealistically large. 

At variance to the previous cases, PBEO and B3LYP 
return now an almost identical adiabatic energy gaps. 
In fact, the B3LYP calculated A£ adia is about 30 meV 
smaller that the PBEO one and, therefore, the trend ob- 
served through the series in Eq. ^ is not respected. As 
we will discuss in detail in the following sections, this 
result might be related to the fact that the energetic of 
[Fe(CO)g] 2+ depends largely the correlation part of the 
functionals as well as the exchange part. Finally we ob- 
serve that the Becke-HH functional incorrectly predict a 
HS ground state, meaning that this includes a too large 
fraction of HF exchange to account accurately for the 
electronic structure of this ion. 



Zero point phononic energies 

So far we have focused only on the adiabatic energy 
gaps. However the expression for the internal energy 
difference, Eq. ([2]), contains also a contribution coming 
from the phonon zero point energies. Table |IX| displays 
AE ZPE , calculated by using the various functionals (the 
results are shown only for the basis set C). Ai? ZPE is 



TABLE IX: Energy difference between the phonon zero point 
energy of the HS and LS state calculated with the various 
functionals employed in this work (only results for the basis 
set C are shown). 



found to be always negative (i.e. the zero point energy of 
the HS state is lower than that of the LS one) reflecting 
the weaker Fe-ligand bond of the HS configuration. Cor- 
rections to the total energy of the two states then always 
tend to stablilise the HS. 

In contrast to the adiabatic energy gap, AE ZPE is 
found to be almost functional independent. Indeed, for 
a given system, the difference between two values of 
AE ZPE obtained with two different functionals, are never 
larger than 15 meV. This demonstrates that the curva- 
ture of the PESs is usually very well reproduced by every 
functional. Therefore the spread in the predicted values 
of A_E adla must arise from the relative shift of the PES 
of one spin state with respect to that of the other. This 
observation is consistent with the results by Zein et al. 
|16) . which indicates that the DOGs, denned by Eq. (p)j), 
do not depend on the choice of functional. 



DMC RESULTS AND DISCUSSION 



The tables |Xl pOj [XTj] and \MU\ display the DMC total 



energies [67] of the four ions [Fe(H 2 0) 6 ] 2+ , [Fe(NH 3 ) 
[Fe(NCH) 6 ] 2 + and [Fe(CO) 6 ] 2+ in both the HS and LS 
states. The molecular geometries were obtained by DFT 
optimization (both the functional and the basis set used 
are indicated in the first column on the left-hand side). 
Unfortunately in most cases the DMC energies have 
a statistical error not small enough to firmly establish 
which functional returns the lowest energy structure of a 
given complex (only the LDA molecular structures have 
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Details geom. opt. 


At (a.u.) 


£ LS (eV) 


-E HS (eV) 




BP86 (Basis C) 


0.0125 


-6127.211(9) 


-6129.720(8) 


-2.51(1) 


BP86 (Basis C) 


0.005 


-6127.218(9) 


-6129.90(2) 


-2.65(1) 


BP86 (Basis C) 


0.001 


-6127.54(9) 


-6130.19(4) 


-2.65(9) 


B3LYP (Basis C) 


0.0125 


-6127.09(2) 


-6129.74(1) 


-2.65(2) 


B3LYP (Basis C) 


0.005 


-6127.36(1) 


-6129.89(2) 


-2.54(1) 


B3LYP (Basis C) 


0.002 


-6127.44(3) 


-6130.01(2) 


-2.57(4) 


B3LYP (Basis C) 


0.001 


-6127.5(1) 


-6130.10(2) 


-2.6(1) 


PBEO (Basis C) 


0.125 


-6127.220(9) 


-6129.804(8) 


2.58(1) 


PBEO (Basis C) 


0.005 


-6127.44(2) 


-6129.94(2) 


-2.50(3) 


PBEO (Basis C) 


0.001 


-6127.66(6) 


-6130.18(4) 


-2.52(7) 



TABLE X: DMC total energy for the LS state, the HS state 
and the adiabatic energy gap of the Fe(H20)e] 2+ ion. The 
molecular structures were optimized by DFT using the vari- 
ous functionals and basis sets listed in the first column. The 
time-steps chosen for the DMC simulation are also indicated. 
Differences in energy are well converged for At = 0.005 a.u. 



Details geom. opt. 


At (a.u.) 


^Ls(eV) 


E HS (eV) 


A£ adia (eV) 


LDA (Basis C) 


0.0125 


-5234.92(1) 


-5236.93(1) 


-2.01(1) 


LDA (Basis C) 


0.005 


-5235.33(2) 


-5237.17(1) 


-1.84(2) 


LDA (Basis C) 


0.001 


-5235.69(5) 


-5237.36(5) 


-1.67(7) 


BP86 (Basis C) 


0.0125 


-5235.56(1) 


-5237.162(9) 


-1.60(1) 


BP86 (Basis C) 


0.005 


-5235.78(1) 


-5237.37(1) 


-1.58(1) 


BP86 (Basis C) 


0.001 


-5235.98(3) 


-5237.55(5) 


-1.57(5) 


B3LYP (Basis C) 


0.0125 


-5235.516(9) 


-5237.15(1) 


-1.63(1) 


B3LYP (Basis C) 


0.005 


-5235.77(1) 


-5237.36(1) 


-1.59(1) 


B3LYP (Basis C) 


0.001 


-5236.01(3) 


-5237.59(4) 


-1.58(5) 


PBEO (Basis B) 


0.0125 


-5235.60(1) 


-5237.21(1) 


-1.61(1) 


PBEO (Basis B) 


0.005 


-5235.89(2) 


-5237.40(2) 


-1.51(2) 


PBEO (Basis B) 


0.001 


-5236.14(3) 


-5237.67(9) 


-1.53(9) 


PBEO (Basis C) 


0.0125 


-5235.57(1) 


-5237.133(8) 


-1.56(1) 


PBEO (Basis C) 


0.005 


-5235.88(2) 


-5237.37(1) 


-1.49(1) 


PBEO (Basis C) 


0.001 


-5236.10(3) 


-5237.60(2) 


-1.50(4) 



TABLE XI: DMC total energy for the LS state, the HS state 
and the adiabatic energy gap of the [Fe(NIi~3)6] 2+ ion. The 
molecular structures were optimized by DFT using the vari- 
ous functionals and basis sets listed in the first column. The 
time-steps chosen for the DMC simulation are also indicated. 
Differences in energy are well converged for At = 0.005 a.u. 



Details geom. opt. At (a.u) 


B LS (eV) 


-EHs(eV) 


A£ ;adia ^yj 


BP86 (Basis C) 


0.0125 


-5957.57(1) 


-5959.30(1) 


-1.73(2) 


BP86 (Basis C) 


0.005 


-5957.57(2) 


-5959.32(2) 


-1.75(3) 


B3LYP (Basis C) 


0.0125 


-5957.94(1) 


-5959.32(1) 


-1.38(2) 


B3LYP (Basis C) 


0.005 


-5957.96(2) 


-5959.33(2) 


-1.37(3) 


PBEO (Basis C) 


0.0125 


-5957.94(1) 


-5959.291(9) 


-1.35(1) 


PBEO (Basis C) 


0.005 


-5957.95(3) 


-5959.30(1) 


-1.35(3) 



TABLE XII: DMC total energy for the LS state, the HS state 
and the adiabatic energy gap of the [Fe(NCH)6] 2+ ion. The 
molecular structures were optimized by DFT using the vari- 
ous functionals and basis sets listed in the first column. The 
time-steps chosen for the DMC simulation are also indicated. 
Differences in energy are well converged for At = 0.0125 a.u. 



systematically higher energies, but this is not surpris- 
ing since the analysis of the phonon modes revealed that 



Details DFT geom. opt. 


At (a.u) 


£ LS (eV) 


£ HS (eV) 


A£ mdia ( e y) 


B3LYP (Basis C) 


0.0125 


-6850.97(2) 


-6850.64(2) 


0.33(3) 


B3LYP (Basis C) 


0.005 


-6850.82(2) 


-6850.45(2) 


0.37(3) 



TABLE XIII: DMC total energy for the LS state, the HS 
state and the adiabatic energy gap of the [Fe(CO)s] 2+ ion. 
The molecular structures were optimized by DFT using the 
functionals and the basis sets listed in the first column. The 
time-steps chosen for the DMC simulation are also indicated. 
Differences in energy are well converged for At — 0.0125 a.u. 



these structures are not even associated to a stable min- 
imum of the LDA total energy) . 



<3 




FIG. 4: Adiabatic energy gaps calculated with DFT and 
DMC. The DFT results were obtained with the functionals 
indicated in the legend and the basis set C. The DMC results 
were obtained for the structures optimized with B3LYP (basis 
set C) and with time-steps At = 0.005 a.u. (the error bars 
are smaller than the symbols). 

In contrast, the adiabatic energy gaps are calculated 
with great confidence and they are listed in the right- 
most column of Tabs. [X XI XII| and |XIII| An analysis 
of these results can be carried out by looking at Fig. |1J 
where we present Ai? adla calculated with both DFT and 
DMC for all the four ions. The systematic up-shift of 
the LDA and BP86 values with respect to the DMC ones 
reflects the massive artificial over-stabilization of the LS 
state (this shift can be as large as few eV) . Notably, LDA 
and BP86 incorrectly return a LS ground state for the 
ions [Fe(NH 3 ) 6 ] 2 + and [Fe(NCH) 6 ] 2 +. 

B3LYP and PBEO provide slightly improved results. 
Their values for AE adla lie systematically below those 
computed with BP86 and the ground state spin is cor- 
rectly predicted for all ions. However, unfortunately, the 
quantitative agreement with DMC is still far from being 
reached as the PBEO and the DMC results differ by about 
0.6 eV (in the best case). Nevertheless hybrid function- 
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3 1 1 1 1 1 1 1 1 1 1 1 1 1— 

2- [Fe(NCH) s r #-#B3LY 

> ■-■ PBEO 

& 1- 




w _ iL [Fe(CO) 6 ] 2+ - 

10 15 20 25 30 35 40 45 
% HF exchange 

FIG. 5: Adiabatic energy gaps versus the fraction of HF ex- 
change included in the hybrid functionals B3LYP and PBEO 
for [Fe(NCH) 6 ] 2+ (upper pannel) and [Fe(CO) 6 ] 2+ (lower 
pannel). The basis set C was used. 



als calculate correctly the relative ligand-field strength 
of the three HS ions [Fe(H 2 0) 6 ] 2+ , [Fe(NH 3 ) 6 ] 2 + and 
[Fe(NCH) 6 ] 2 +. In fact, although the B3LYP and PBEO 
results appear shifted vertically according to the fraction 
of HF exchange included in the functional, the relative 
Ai5 adla of two complexes is well reproduced. In contrast, 
the results for [Fe(CO)g] 2+ do not show the same trend. 
The PBEO adiabatic energy gap lie slightly above (by 
about 30 meV) the B3LYP one, despite the fraction of 
HF exchange being larger in PBEO than in B3LYP. This 
indicates that the exchange and correlation energies have 
different relative importance for the HS and the LS com- 
pounds. In order to better understand this important ob- 
servation, we have calculated the adiabatic energy gap of 
[Fe(NCH) 6 ] 2 + and [Fe(CO) 6 ] 2+ after changing the frac- 
tion of HF exchange in both B3LYP and PBEO (see Fig. 

§■ 

On the one hand, PBEO and B3LYP give very simi- 
lar results for [Fe(NCH)g] 2+ , regardless of the different 
local-exchange and correlation energy. Therefore the cal- 
culated adiabatic energy gaps depend mainly on the frac- 
tion of HF exchange (and this dependence is almost lin- 
ear). This indicates that the correlation contribution to 
the total energy is well described with (semi-)local func- 
tionals and the failures in predicting A_E adla could be 
entirely ascribed to the underestimation of the exchange 
energy. In addition, by fitting the data, we also conclude 
that about 50% of HF exchange is required to achieve a 
fair agreement between the DFT and the DMC adiabatic 
energy gaps. Hence, the Becke-HH functional is found to 
provide quite satisfactory results (see Fig. [I]). 

On the other hand, the B3LYP-calculated A£ adia of 
[Fe(CO) 6 ] 2+ is systematically down-shifted with respect 
to the PBEO value calculated with the same amount of 
HF exchange. Therefore for this ion the results depend 



drastically on the correlation as well as the exchange part 
of the density functional. In addition, we note that about 
30% and 40% of HF exchange, respectively in B3LYP and 
PBEO, is required to reproduce the DMC gaps and that 
the HH functional incorrectly describes [Fe(CO)e] 2+ as a 
HS ion (see Fig. [I]). 

The DFT performances for the three HS ions and 
[Fe(CO) 6 ] 2+ is related to the different nature of their 
ground state wave-function. In fact, this was found 
to have a much more pronounced multi-configurational 
character in [Fe(CO)e] 2+ than in all the other complexes 
[68] , reflecting the increase in the covalency of the metal- 
ligand bonds [551 IBS]- Based on this observation, one can 
reasonably argue that, for the HS complexes the local- 
part of the exchange-correlation functional is able to cap- 
ture most of the correlation energy, while the fraction of 
HF exchange effectively cures the LDA underestimation 
of the exchange. Thus hybrid functionals with "enough" 
HF exchange are found to systematically return quite 
satisfactory results. Furthermore, as the failures of stan- 
dard GGA functionals seem mostly related to the short- 
comings in the description of the exchange energy, recent 
GGA functionals, constructed in order to tackle this is- 
sue, can out-perform. For example, the OLYP functional, 
whose exchange part (OPTX) is parametrized to repro- 
duce the Hartree-Fock exchange for atoms [23], predicts 
values for the adiabatic energy gap of [Fe(H 2 0)6] 2+ and 
[Fe(NH3)g] 2+ , which compare well with those computed 
with cither B3LYP or PBEO gS]. Unfortuantelly, how- 
ever, OLYP is not as accurate as the hybrids for predict- 
ing geometry optimizations and bond-lengths [46 . An- 
other of such GGA functionals, which was found to per- 
form at the level of the hybrids (if not even better) [26] , 
is the HCTH407 [69]. Our own results are then listed in 
Tab. |XIV| The massive improvement, that OLYP and 
HCTH407 have achieved, over BP86, is evident in the 
case of Fe(H 2 0) 6 ] 2+ , [Fe(NH 3 ) 6 ] 2 + and [Fe(NCH) 6 ] 2 +. 



System 


Functional Basis set AE 


idia (cm" 


-i) A£ adia (eV) 


[Fe(H 2 0) 6 ] 2+ 


OLYP 


C 


-15953 


-1.9780 


[Fe(H 2 0) 6 ] 2+ 


HCTH407 


C 


-19315 


-2.3947 


[Fe(NH 3 ) 6 ] 2+ 


OLYP 


C 


-7338 


-0.9099 


[Fe(NH 3 ) 6 ] 2+ 


HCTH407 


c 


-9942 


-1.2327 


[Fe(NCH) 6 ] 2+ 


OLYP 


c 


525 


0.06510 


[Fe(NCH) 6 ] 2+ 


HCTH407 


c 


-3650 


-0.4526 


[Fe(CO) G ] 2+ 


OLYP 


c 


21313 


2.6425 


[Fe(CO) 6 ] 2+ 


HCTH407 


c 


17097 


2.1198 



TABLE XIV: Adiabatic energy gap, A_E adia , for the four ions 
calculated with the OLYP and HCTH407 functionals (the ba- 
sis sets C was used). 

In contrast, one can question whether [Fe(CO)g] 2+ can 
be described at all by the single-determinant picture pro- 
vided by DFT. In principle, the multiconfigurational na- 
ture of a wave-function can be described by GGA func- 
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tionals. In fact the GGA exchange roughly mimics the 
non-dynamical correlation (in addition to the proper ex- 
change) [70]. In practice, however, no DFT flavour inves- 
tigated here has proven fairly accurate for the energetic 
of the ion [Fe(CO) 6 ] 



12+ 



Method 


Reference 


A£ adia (cm" 1 


) A£ adia (eV) 


CASSCF(6,5) 


[25] 


-23125 


-2.86714 


CASSCF(12,10) 


1251 


-21180 


-2.62599 


corr-CASSCF(12,10) 


m 


-17892 


-2.21833 


CASPT2(6,5) 


E3 


-21610 


-2.6793 


CASPT2(12,10) 


EU 


-16185 


-2.00668 


corr-CASPT2(12,10) 


.251 


-12347 


-1.53083 


SORCI 


251 


-13360 


-1.65643 


CASPT2(10,12) 


[651 


-16307 


-2.02181 



TABLE XV: The adiabatic energy gap for [Fe(H 2 0) 6 ] 2+ cal- 
culated by using various wave-function methods (reference to 
the literature is given in the second column). The values la- 
belled as corr-CASSCF and corr-CASPT denote respectively 
the CASSCF and CASPT values after having applied an em- 
pirical correction of the order of 4000cm -1 (see main text). 
Pierloot et al. [65| provides an additional long list of results 
obtained by using different basis sets, geometries and sym- 
metries. Here we report only the value that these authors 
indicate as the "best" . 



These were stated not to require any empirical correc- 
tion. A second set of results is provided by Pierloot et al. 
[65] , who performed calculations with basis sets of larger 
size. For [Fe(NH 3 ) 6 ] 2+ , their best CASPT2 adiabatic 
gap agrees fairly well with the corrected-CASPT2 and 
SORCI results. However for [Fe(H 2 0) 6 ] 2+ , they found a 
significantly larger (in absolute value) Ai? adla . 

we note 



By analyzing the data in Tabs. |XV| and XVI 
that the adiabatic energy gaps calculated with CASPT2 
by Fouqueau et al. [35] [25] agree fairly well with our 
DMC ones (in particular for [Fe(NH 3 ) 6 ] 2+ ). In contrast, 
the empirical corrections worsen the agreement and the 
SORCI results do not agree quantitatively with ours. Al- 
though we have not achieved yet a complete understand- 
ing of these differences, we argue that they may originate 
from the large dependence of the CASSCF/CAST2 re- 
sults on the basis sets and on the orbitals included in the 
active space. Furthermore wave-function based methods 
do not describe dynamic electronic correlations (although 
partial corrections are provided by the second order per- 
turbation theory). The DMC energies, in turns, might 
depend on the choice of the trial wave-function intro- 
duced to impose the fixed-node approximation [37]. A 
thorough analysis on the sources of potential errors in 
DMC is currently on the way. 



Method 


Reference 


A£ adia (cm" 1 ) 


A£ adia (eV) 


CASSCF(12,10) 


E<D 


-20630 


-2.55779 


corr-CASSCF(12,10) 


Ea 


-16792 


-2.08194 


CASPT2(12,10) 


261 


-12963 


-1.60721 


corr-CASPT2(12,10) 


261 


-9125 


-1.13136 


SORCI 


Eg - 


-10390/ - 11250 


-1.2882/ - 1.39482 


CASPT2(12,10) 


1551 


-7094 


-0.879544 


TABLE XVI: Adiabatic ener 


gy gaps for [Fe(NH3)e] 2+ calcu- 



lated by using various wave-function methods (reference to 
the literature is given in the second column). The values la- 
belled as corr-CASSCF and corr-CASPT denote respectively 
the CASSCF and CASPT values after having applied the em- 
pirical correction of the order of 4000cm -1 (see main text). 

Finally we compare our DMC results to those obtained 
with wave-function based methods. Fouqueau et al. 
[231 l2"rj] carried out several calculations for the adiabatic 
energy gap of the ions [Fe(H 2 0) 6 ] 2+ and [Fe(NH 3 ) 6 ] 2+ 
by using the complete active space self-consistent field 
(CASSCF) method with second order perturbative cor- 
rections (CASPT2). Some of the results are summarized 
in Tabs . [XV] and |XVl] As observed by the authors them- 
selves and in Ref . [53] , these calculations suffer the draw- 
back of having been carried out with an insufficient Fe 
basis set. As such they are affected by a systematic er- 
ror, which can be estimated by considering the 5 D — 1 / 
splitting of the free Fe 2+ ion. An empirical correction 
of the order of 4000 cm -1 was then introduced. In the 
same works, results obtained by spectroscopy-oriented 
configuration-interaction (SORCI), were also reported. 



CONCLUSIONS 

In this work, we have assessed the performances of sev- 
eral popular exchange-correlation functionals in describ- 
ing various Fe 2+ complexes. As DFT results can not be 
easily related to experiments (at least without account- 
ing for environmental and finite-temperature effects), we 
have performed accurate DMC calculations, which pro- 
vide a solid theoretical benchmark for the theory. The 
DFT and DMC results, both obtained within the theo- 
retical framework of the adiabatic approximation, could 
be then directly compared. 

The LDA and the standard GGA functionals drasti- 
cally over-stabilize the LS state. Although the accuracy 
of the DFT calculations increases when hybrid function- 
als are employed, the most popular ones, B3LYP and 
PBE0, provide results, which are still quantitatively un- 
satisfactory. In the case of HS ions, a fair agreement be- 
tween the DFT and the DMC adiabatic gaps is achieved 
only by using about 50% of HF exchange. In contrast, a 
lower fraction of HF exchange (between 30% and 40%) 
is required for [Fe(CO)g] 2+ . This difference might be 
related to the diverse nature of the ground state wave- 
function for the HS and LS ions. Therefore, unfortu- 
nately, we have to conclude that there is not yet a "uni- 
versal" functional able to correctly describe the energet- 
ics of every Fe 2+ complex. 

Finally, by analyzing zero-point phonon energies, we 
have demonstrated that the shape of the PESs is well 



12 



described by every functional considered. Therefore, as 
already pointed out by Zein et al. [16], the failures of 
DFT in calculating the adiabatic energy gaps must be 
ascribed to a shift of the PES of the LS state with respect 
to that of the HS state. 
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